Lattice Boltzmann model for ultra-relativistic flows 
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We develop a relativistic lattice Boltzmann model capable of describing relativistic fluid dynamics 
at ultra-high velocities, with Lorentz factors up to 7 ~ 10. To this purpose, we first build a new 
lattice kinetic scheme by expanding the Maxwell- J iittner distribution function in an orthogonal 
basis of polynomials and applying an appropriate quadrature, providing the discrete versions of the 
relativistic Boltzmann equation and the equilibrium distribution. To achieve ultra-high velocities, we 
include a flux limiter scheme, and introduce the bulk viscosity by a suitable extension of the discrete 
relativistic Boltzmann equation. The model is validated by performing simulations of shock waves 
in viscous quark-gluon plasmas and comparing with existing models, finding very good agreement. 
To the best of our knowledge, we for the first time successfully simulate viscous shock waves in 
the highly relativistic regime. Moreover, we show that our model can also be used for near-inviscid 
flows even at very high velocities. Finally, as an astrophysical application, we simulate a relativistic 
shock wave, generated by, say, a supernova explosion, colliding with a massive interstellar cloud, 
e.g. molecular gas. 

PACS numbers: 47.11.-j, 02.40.-k, 95.30. Sf 



I. INTRODUCTION 

Relativistic fluid dynamics plays an important role in 
many contexts of astrophysics and high-energy physics, 
e.g. jets emerging from the core of galactic nuclei or 
gamma-ray bursts pQ , shock induced Ritchmyer-Meshkov 
instabilities [2] and quark-gluon plasmas produced in 
heavy- ion collisions [3]. Hence, various numerical meth- 
ods have been developed to study the relativistic hy- 
drodynamics. Most of these methods are focussed on 
the solution of the corresponding relativistic macroscopic 
conservation equations. Among others, one can men- 
tion the methods based on second-order Lax-Wendroff 
scheme [J, smoothed particle hydrodynamics techniques 
[5j [6] , Glimms (random choice) method [7] and high res- 
olution shock-capturing methods [8]. Other methods, in- 
stead of solving the macroscopic equations, tackle the 
problem from the microscopic and mesoscopic points of 
view [9 j. To this regard, the lattice Boltzmann (LB) 
method jT0UT2] a relatively new numerical approach, 
based on a minimal lattice version of the Boltzmann ki- 
netic equation, has enjoyed increasing popularity for the 
last two decades. Within LB, representative particles 
stream and collide on the nodes of a regular lattice, with 
sufficient symmetry to reproduce the correct equations of 
macroscopic hydrodynamics. The main highlights of LB 
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are its computational simplicity, easy handling of com- 
plex geometries, and high amenability to parallel com- 
puting [13]. The LB method has met with remarkable 
success for the simulation of a broad variety of complex 
flows, from fully developed turbulence, all the way down 
to nanoscale flows of biological interest [T4HT6] . 

From a mathematical viewpoint, the standard lattice 
Boltzmann model can be obtained by expanding the equi- 
librium distribution, i.e. Maxwell-Boltzmann distribu- 
tion, in a Hermite polynomials and using the nodes of 
polynomials, up to a certain order as the correspond- 
ing discretized velocities [17], using the Bhatnagar-Gross- 
Krook (BGK) approximation for the collision operator 
[T8] . While the applications of the LB scheme cover an 
impressive array of complex fluid flows, its relativistic ex- 
tension has been developed only in last few years p~9j [20] . 

The relativistic LB (RLB) model was constructed by 
expanding the distribution function in powers of the 
fluid speed and finding the corresponding coefficients 
(Lagrange multipliers), by matching the moments of 
the Maxwell- Jiittner distribution in continuum velocity 
space. This model was shown capable of simulating 
weakly and moderately relativistic viscous flows, with 
/3 = ii/c~0.3, u being the typical flow speed. In partic- 
ular, RLB was applied to the simulation of shock waves 
in quark-gluon plasmas, showing very good agreement 
with the results obtained by solving the full Boltzmann 
equation for multi-parton scattering (BAMPS), [21], 

However, the aforementioned matching procedure does 
not provide a unique solution for the discrete equilib- 
rium distribution function, satisfying the hydrodynam- 
ics moments of the Maxwell- Jiittner distribution. More- 
over, the model lacks dissipation for the zero component 
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of the energy-momentum tensor, and imposes a non- 
physical diffusion in the conservation of the number of 
particles [22 . These flaws, albeit very minor at moder- 
ate flow speeds, may become a concern for strongly rel- 
ativist ic flows. It is therefore highly desirable to develop 
more general and systematic approaches. To this pur- 
pose, let us observe that, due to the non-separability of 
the Maxwell- Jiittner distribution function into the three 
components of the momentum in Cartesian coordinates, 
its expansion in orthogonal polynomials is not as natu- 
ral as in the classical case and some deliberation is re- 
quired. For the fully relativistic regime, neglecting par- 
ticle masses, and by using spherical coordinates, a lat- 
tice Boltzmann algorithm for the relativistic Boltzmann 
equation was developed in Ref. [23]. In this paper, the 
Maxwell- Jiittner distribution function was expanded in 
an orthogonal polynomials basis and discretized using a 
Gauss quadrature procedure. The model was based on 
the Anderson- Witting collision operator [24] . The results 
of simulating viscous quark-gluon plasma were compared 
to other hydrodynamic simulations and very good agree- 
ment was observed. However, using spherical coordinates 
makes the scheme incompatible with a cartesian lattice, 
and consequently, in the streaming procedure, a linear 
interpolation is required at each time step. Therefore, 
some crucial properties of the classical LB, e.g. exact 
streaming (zero numerical dispersion) and negative nu- 
merical diffusivity, are lost in the process. 



II. MODEL DESCRIPTION 

We start the description of our model by writing the 
Maxwell- Jiittner equilibrium distribution function as 



n = AeM~P^/k B T), 



(1) 



where, A is a normalization constant, k B the Boltzmann 
constant, T the temperature and = (E/c,p) is the 
4-momentum, with the energy of the particles E defined 

by 



E = cp° 



mc 



a/1 - u 2 /c 2 



(2) 



The macroscopic 4-velocity is (U^) = (c,u)j(u), with u 
the three-dimensional velocity, j(u) = 1/yT — u 2 /c 2 the 
Lorentz's factor, m the mass, and c the speed of light. 
The relativistic Boltzmann equation, based on the Marie 
collision operator [25 , reads as follows 



: (/"/ eq ), 



(3) 



where / is the probability distribution function, and r the 
single relaxation time. It is possible to write the Maxwell- 
Jiittner distribution in a simpler form by introducing the 
following change of variables: 



(4) 



In this paper, we develop a relativistic lattice Boltz- 
mann model by expanding the Maxwell- Jiittner distribu- 
tion in a set of orthogonal polynomials, and performing 
an appropriate quadrature in order to adjust the scheme 
to a D3Q19 (19 discrete velocities in three spatial dimen- 
sions) cell configuration [19 , 20 . Moreover, we extend the 
model by using a minimum modulus flux limiter scheme 
and introducing the bulk viscosity term into the Boltz- 
mann equation. We show that the model is numerically 
stable also at very high velocities, i.e. Lorentz factors up 
to 7 ~ 10. Additionally, we show that this model can 
also be used to simulate near-inviscid flows, which cor- 
responds to solve the Euler equation on the macroscopic 
level. This is well suited for astrophysical applications, 
where the viscosity is usually negligible. In fact, the as- 
trophysical context presents possibly the richest arena for 
future applications of the present RLB scheme. 



The paper is organized as follows: in Sec. pl} th e model 
description is presented in detail; and in Sec. |lll[ several 
validations with other existing numerical models along 
with some results for shock waves in viscous quark-gluon 
plasmas and a 3D simulation of a shock wave colliding 
with a massive interstellar cloud are presented. Finally, 
in Sec. |IV[ a discussion about the model and the results 
is provided. 
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m 



v = c/c s , 



(5) 



and therefore, by replacing the new variables in Eq. ([!]) 
we have 



r = AeM-^xn- 



(6) 



The temporal components, £° and x°? can be calculated 
by the relations 



(7) 



(8) 



In analogy to the classical procedure of expanding the 
Maxwell-Boltzmann distribution in Hermite polynomi- 
als, we can also expand the Maxwell- Jiittner distribution, 
using orthogonal polynomials of the following form: 
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where 
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for m^n, and 
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N {n) = I wFi^F, 
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(n) r {n) £ • 



(11) 



To construct the appropriate orthogonal polynomials, we 
introduce the corresponding weight function as the equi- 
librium distribution at the local rest frame, 



w(£) = Aexpi-vt ). 



(12) 



Using the procedure proposed by Stewart [26], where 
the non-equilibrium distribution was expanded around 
the equilibrium, and the Maxwell- Jiittner distribution 
was used as the weight function to find the orthogonal 
polynomials, we can take up to second order, 



and 



F (o) = 



I'm — C - a , 



(13) 



(14) 



pap = ^ - a fF^ - b a ?, (15) 

where a", a a ^ and b a ^ are unknowns to be calculated 
using the Gram-Schmidt orthogonalization procedure 



/ 



The normalization coefficient for each polynomial is given 
by y/7V( n ) , and the coefficient is calculated using the 
relation 



a (n) = I / eq ^(n) 



(17) 



To calculate the coefficients a a , a a/3 and 6 a/3 , one needs 
the moments of the Maxwell- Jiittner distribution, which 
up to second order are given by [27] 



J n d ^=A^AK 1 (v% 



(18) 
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r/ eq = ^AK 2 (u 2 ) X a , (19) 
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r^/ eq § = {K 2 { V 2 )^ - K 3 (u 2 ) X a X ) , 



£0 



(20) 



where K n (y 2 ) is the modified Bessel function of the sec- 
ond kind of order n and r] a ^ is the Minkowski metric 
tensor with the signature (+,—,—,—). The moments 
with respect to the weight function can be determined 



by considering the above integrals in the local Lorentz 
rest frame. 

For the sake of simplicity we define <p a as 



l ) = (x°,o), 



(21) 



and by using the orthogonalization relations to calculate 
the unknowns, the resulting relativistic orthogonal poly- 
nomials are given by 



(o) 



1, 



rpa tot 



rpOtfi 
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where 



h aP 



K 3 (y 2 ) 
K^p 2 )' 
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K 2 (y 2 ) 
K^p 2 ) 



(22) 
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and 
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- (k 3 (u 2 ) 

+ (k a {v 2 
Here, the function D{y) is defined by 

[i + ■ K ^ 2) 



(26) 
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By following this procedure, we can calculate higher 
order polynomials. However, since in this work we are 
interested in recovering only up to the second moment 
of the Maxwell- Jiittner distribution (energy-momentum 
tensor), using the expansion up to the second order is 
sufficient. In particular, the third, fourth and fifth or- 
der moments, which are needed to describe highly vis- 
cous fluids, would increase dramatically the complexity 
of our expansion, and consequently its numerical imple- 
mentation. This is a very interesting subject for future 
extensions of this works. 

Additionally, in the ultrarelativistic limit, where 
ksT ^> mc 2 , i.e. v <C 1, we can use the following asymp- 
totic relation: 

2 n - 1 (n- 1)! 



(28) 



Using the resulting polynomials -F( n )> coefficients a( n ) 
and iV( n ) with Eqs. (|9|) and (28) we can expand the 



Maxwell- Jiittner distribution in orthogonal polynomials, 

H x Z v x x x v + t x t z x x x z + t v ¥x v x z 

— 4 — [(z x )Hx x ) 2 + m 2 {x v ? + (e) 2 (x z ) 2 



(29) 
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FIG. 1. Comparison between the Maxwell Juttner distribu- 
tion function and the zeroth, first and second order expansions 
in one space dimension for f3 = 0.2. In the inset, i 7 ^) 

ploynomials, in the local rest frame are 
shown. 



Ffa and Fffi 



up to second order and in the ultrarelativistic regime. 

One can compare the Maxwell Juttner distribution 
with the zeroth, first and second order expansions in the 
one dimensional case. The result of the distributions ver- 
sus £ x for the case j3 = \u\/c = 0.2 is presented in Fig. 
[1] Here, we can observe that, as expected, as the order 
of the expansion increases, the expansion becomes more 
accurate. Note that the expanded distributions become 
negative for values of £ x , around —3. However, this is of 
no concern for our model since, as we shall see shortly, 
the quadrature requires only £ x ~ 1 . For illustrative pur- 
poses, in the inset of Fig. [I] we show, for v = 1 and in 
the local rest frame, the polynomials corresponding to 
the zeroth (-F( )), first (F^), and second (F^j) orders. 

We can write the Boltzmann equation, Eq.Q, as fol- 
lows: 



rc 



r 



(30) 



where latin subscript a runs over the spatial coordinates. 
In order to discretize Eq.(30) and avoid a multi-time lat- 



tice, we need to consider the temporal components of the 
discretized velocity 4- vector, i.e. <^°, as constant. There- 
fore, a transformation of the temporal component of both 
£ a and x a i s required. We can write (£f) = (c t /co,c a ) 
and (x a ) = (x°/ c o,x)i where q, cq and c a are constants 
related to the size of the lattice. We will use a lattice 
configuration D3Q19 (19 discrete vectors in 3 spatial di- 
mensions), which can be expressed as 



(0,0,0) i = 0; 

c a (±l,0,0) F5 1<^<6; 
c a (±l,±l,0) F5 7<z<18, 



(31) 



where the subscript FS denotes a fully symmetric set of 
points. 



To find the discretized weights for the lattice, we use a 
quadrature procedure. According to the quadrature rule, 
the discretized weights should satisfy the relation 



/ 



N 



(32) 



where R(£) is an arbitrary polynomial of order 27V or 
less. Using this relation, we can construct a system of 
equations by replacing R(£) with different combinations 
of zeroth, first and second order polynomials. The left 
hand side of the above equation can be calculated by us- 
ing Eq.(18) to (20). Thus, the resulting discrete weights 



are given by 



w = l 



4c> 2 



2 n 2 ' 



2166c 2 c 2 a 



(361 



a°0 



8cV) 



for 1 < i < 6, and 



Wi = 



1083^ ' 



(33) 



(34) 



(35) 



for 7 < i < 18. Note that we still need to calculate the 
constants related to the size of the lattice, i.e. c ai c t and 

Co- 

The discretized 4-momenta should satisfy the following 
relation for the energy-momentum tensor, i.e. the second 
order moment of the distribution function, 
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d 3 P 



N 



p 1=1 
= (e + p) 
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prj 



a(3 



TSt , (36) 



where p is the hydrostatic pressure, e the energy density 
and T%fi denotes the energy-momentum tensor at equi- 
librium. Note that higher order moments of the discrete 
equilibrium distribution can be calculated by perform- 
ing the respective sums T a/3 - - 7 = Y^iLi P?Pi —Pift^- 
However, they would not correspond to the ones of the 
Maxwell- Juttner distribution, because the latter require 
an expansion in higher order polynomials. Their contri- 
bution to the dynamics of the fluid become important 
at high values of the Knudsen number (high momentum 
diffusivity) , and since they are not exactly recovered, our 
model does not work properly in that regime. Fortu- 
nately, many applications in astrophysics and high energy 
physics deal with near-inviscid or weakly viscous fluids. 

We can simply find the constants related to the lattice 
size, using the fact that in the tensor, the coefficient of 
jjajj(3 £ Qr different a and j3 should be always the same. 
The calculated values for the constants are 



19 



ct/co 



27 



\(9-2y/3). (37) 
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In the ultrarelativistic limit and considering the natu- 
ral units c = ks = 1, from the energy- momentum tensor, 
one can obtain the following relations: 



P 



4n 
7^2 ' 



P 



n 

9 ' 



3p 



(38) 



finding that the relation between e and p corresponds 
to the well-known state equation in the ultrarelativistic 
limit. 

Note that due to the fact that we have supposed 
to be constant, to avoid a multi-time evolution lattice, 
there are some equations in the quadrature procedure 
for the first order moment and the second order moment 
of the distribution function which could not be satisfied 
simultaneously. Indeed, we can choose whether to re- 
cover the first order moment or the second order moment 
in the quadrature. To satisfy the first moment of the 
Maxwell- Jiittner distribution function leads to recover 
the equation for the conservation of number of particles, 
d a N a = 0, and the second order moment, the equation 
for the conservation of momentum-energy, d a T a P = 0. 
To calculate the four unknowns, namely U x ,U y , U z and 
e, the four equations corresponding to T 00 ,T 0x ,T 0y and 
T 0z components of the energy-momentum tensor and 
equation of state e = 3p would be enough. Therefore, by 
using the ultrarelativistic equation of state the dynamics 
of the system is not affected by the number of parti- 
cles and the equation for the conservation of momentum- 
energy is therefore sufficient to describe the entire dy- 
namics of the relativistic fluid. For this reason, our 
quadrature is targeted to recover the second order mo- 
ment, using a separate distribution function to recover 
the equation for the conservation of number of particles, 
d a N a = 0, based on the model proposed by Hupp et al. 



At each time step and at each lattice point, the val- 
ues of the macroscopic velocity and energy density can 
be evaluated using the energy-momentum tensor as men- 
tioned previously. 



A. Extended model for high velocities 

At high velocities (/3 > 0.6), due to the compressibil- 
ity effects (high Mach numbers), the described numeri- 
cal scheme shows artificial discontinuities in the velocity 
and pressure profiles, leading to numerical instabilities in 
the long-term evolution. We shall return to this issue in 



Sec. Ill The relativistic Mach number can be expressed 
as M K — j(u)\u\/j(c so )c so where c so is the velocity 
of sound, which is c so — cf y/3 in the ultra-relativist ic 
regime. In order to overcome this problem, we first use a 
modified version of the D3Q19 cell configuration, which 
is denoted by = (ct/co, c a ), where c a takes now the 
following form: 



(0,0,0) z = 0; 

c a (±l,0,0) FS 1 < i < 6; 
2c a (±l,±l,0) F5 7 < z < 18. 



(42) 



Since some of the discrete velocities go beyond first and 
second neighbors in the lattice, the scheme can sup- 
port higher flow speeds. In order to find the discretized 
weights, as well as the size of the lattice cells, we repeat 
the same procedure as before obtaining 



w = 1 + 



7cj^ 
676cg 



2 r 2 ' 



(43) 



a°0 



Using the mentioned lattice to discretize the Boltz- 
mann equation, Eq.(30) can be written as follows: 



Co 



fi(x + Ca-±5t,t + 5t)-fi(x,t) 



C^vSt 

rc t 



(fi-fr). 



(39) 

The left hand side of the equation is readily recognized as 
free-streaming, while the right hand side is the discrete 
version of the collision operator according to the model 
of Marie. In this equation, the following relation between 
St and Sx has been used: 



St-- 



c t Sx 
c a co ' 



(40) 



In the ultra-relativistic limit, the shear viscosity using 
the model of Marie can be calculated as: 



V : 



(e + p) 



(41) 



In our model, this expression is only valid for small 
values of r (which leads to small values of the Knudsen 
number), where higher order moments of the distribution 
can be neglected. 



Wj = 



2 v 2 ) 



for 1 < i < 6, 



w 1 



8112c§' 



for 7 < i < 18, and 



C a 



(44) 



(45) 



— , c £ /c = ^, c = l(9-2V3). (46) 
v v o 



The second step to extend our model, is to introduce 
the minimum modulus (min mod) scheme to discretize 
the spatial components in the streaming term in the 
Boltzmann equation, Eq.(30), i.e. pfd a fi. The min mod 



scheme is a flux limiter method that efficiently reduces 
the instability especially when step discontinuities occur 
(e.g. in shock waves). The following relations character- 
ize this scheme [29l: 



daiptfi) = 



1 



\Sxe a \ 



[h a i {x + 8xe a )-h a i {x)]. 



(47) 
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-^min mod 



K{x) = ft (*) + /? (x), 



ff{x) = ff{x) 
(^Aff(x),Aff(x-5xe a 



ff(x) = ff(x + 5xe a ) 
- imin mod ( A/f + (x) , Aff (x + 5xe a )) , 



(48) 



(49) 



(50) 



ff = \{p1 + \pt\)h, ff = \{p a i-\p a i\)fi, (si) 



A/f ± (x) = ff (x + Sxe a ) - ff (x), 



(52) 



where e a is a unit vector in the direction of the corre- 
sponding spatial coordinate. Let us remind that is 
independent of spatial coordinates. The min mod func- 
tion is defined as 



min mod(X,r) = |min(|X|, \Y\) 
x[Sign(X) + Sign(y)]. 



(53) 



Note that in classical flows, the bulk viscosity plays 
an important role in highly compressible flows and en- 
hances the stability of numerical simulations of fluids at 
very high velocities, including shock waves [30]. How- 
ever, in the ultrarelativistic limit, the energy- momentum 
tensor is traceless and the bulk viscosity is zero [31] [32] . 
Therefore, in order to include the bulk viscosity, we add 
the following term to the right hand side of the Boltz- 
mann equation, Eq.([30|): 



(54) 



where 



a=x,y,z 







i 0; 
a5x i 7^ 0, 



(55) 



where a is a constant. A central finite difference 
scheme is used to calculate the second order derivative. 
Chapman-Enskog analysis reveals that the bulk viscosity 
obtained by this extra-term takes the form 



rib 



4pau 6 
27 



(56) 



Note that, like the analytical expression of the bulk vis- 
cosity for the model of Marie, the above-mentioned bulk 
viscosity is also proportional to T -3 , and as expected 
goes to zero in the ultra-relativistic limit v — »> 0. How- 
ever, this small value of bulk viscosity is sufficient to 
stabilize the system at high velocities, as we are going to 
show in the next section. 



Once all of the extensions above are taken into ac- 
count, the discretized form of the relativist ic Boltzmann 
equation takes the following expression: 



cpvdt 
rc t 



fi(x,t + 5t)-Mx,t) 
■%g[h?(x + 5xe a )-h*(x)] = 



[fi(x,t 



-j*)-(2/rw)-/rw- 

- c -^\iT,dlfi{x,t), 



St))] 



(57) 

where an implicit representation of the collision term is 
used, as proposed in Ref. [33 to enhance the stability of 
the collision term. 

As mentioned above, for the cases where the dynamics 
of the number of particles density is also needed, one 
has to solve the conservation equation, i.e. d a N a = 0, 
with N a = nU a . For this purpose, we add an extra 
distribution function, hi, based on the model proposed 
by Hupp et al. [28], which follows the dynamics of the 
Boltzmann equation given by Eq. (57), without the 



coefficient term. The corresponding modified equilibrium 
distribution function is given by: 



h? = w'n^u) ( | + 3(c a .u) + \(c a .uf - \\u\ 



Wn 



1 

To' 



_6_ 
35 



42cT 



for 1 < i < 6, 



84c2 



3 
280 



(58) 
(59) 



(60) 



for 7 < i < 18, where we have used our new cell configu- 
ration, and w[ are the respective discrete weights. 

Having discussed the model, we move on to the next 
section, where different validations and results for the 
Riemann problem are provided, along with a simulation 
of a shock wave colliding with a interstellar cloud. 



III. VALIDATION AND RESULTS 

In order to validate the model and the numerical proce- 
dure, we present results for the simulation of relativist ic 
shock wave propagation in viscous quark-gluon plasma. 
The associated Riemann problem is studied and several 
comparisons are drawn between the present RLB model 
and the existing literature. Indeed, the Riemann prob- 
lem is a challenging test for numerical methods, since it 
involves a shock and rarefaction wave. 

The initial condition of the Riemann problem consists 
of two regions with different pressure, which are sepa- 
rated by a membrane in the middle of the interval. The 
pressure in the left region (P ) is higher than the pressure 
in the right region (Pi). Both sides of the discontinuity 
are supposed to be initially in the rest frame. Hence, spa- 
tial components of initial velocities for both sides are set 
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FIG. 2. Comparison of the velocity profile of the current 
model and BAMPS, for different values of r]/s at weakly rel- 
ativistic regime. 



to zero. At time t = 0, the membrane is removed and a 
shock wave propagates from the high pressure region into 
the low pressure region with velocity v s hock and a rar- 
efaction wave propagates in the opposite direction. The 
shock velocity only depends on the pressure difference, 
the equation of state, and can be calculated analytically 
[34| [35] . The region between the shock wave and the rar- 
efaction wave has a constant pressure, corresponding to 
the so-called shock plateau. In this region the velocity is 
also constant (v p i at ). 

In order to compare our results with existing models, 
we use the same conditions as in Ref. [HI [21] . There- 
fore, one dimensional simulations are carried out using 
800 x 1 x 1 cells, where open boundary conditions are 
considered at the two ends. The cell size S x is taken to 
be unity, which corresponds to S x = 0.008/m in IS units 
and St can be calculated from Eq.(40). We use n/s as 
the characteristic parameter of shear viscosity, where the 
viscosity is defined in Eq.(41) and the entropy density is 
given by s = An — nlnA, with A = n/n eq being the fu- 
gacity of gluons, n eq = dcT 3 /tt 2 the equilibrium density, 
and do = 16 the degeneracy of gluon. 

For the first validation test, we set the initial pres- 
sure at the left and right sides to Po = 5A3GeV/ fm 3 
and Pi = 2.22GeV/ fm 3 , respectively. This corresponds 
to 2.495 x 10" 7 and 1.023 x 10" 7 in numerical units, 
respectively. Fig. [2] shows the pressure profiles at time 
t = 3.2 fm/s for different values of 77/s, compared to the 
results reported by Ref. [21] (hereafter BAMPS) and an 
analytical solution for the inviscid case reported in Ref. 
[34] . As one can notice, very satisfactory agreement for 
different values of n/s is obtained.. 

As mentioned previously, the lattice Boltzmann 
method is computationally very efficient. For instance, 
the above simulation took ~ 220 ms on a standard PC, 
which is approximately an order of magnitude faster 
than corresponding hydrodynamic simulations. To fur- 
ther elaborate on the validity of our model, we compare 
with the results of Ref. [19] (hereafter previous LBS) for 



Current Model: r|/s=0.005 - 
Previous LBS: r|/s=0.005 
Current Model: r|/s=0.01 

Previous LBS: r|/s=0.01 
Current Model: r|/s=0.05 - 

Previous LBS: r|/s=0.05 




FIG. 3. Comparison between the current model and the pre- 
vious LBS model at different 77/s, for the pressure (top) and 
velocity (bottom) profiles in the weakly relativistic regime. 



different values of n/s. Fig. [3] shows that pressure and ve- 
locity profiles are in good agreement with previous LBS 
simulations. It is worth mentioning that, as it is apparent 
from Fig. [3] the above mentioned pressure difference cor- 
responds to f3 = \vpi at \ ~ 0.2 (weakly relativistic regime), 
while the velocity of the shock is t? s ho C fc ~ 0.65. 

To study higher velocities, we consider higher pres- 
sure difference between the left side and the right side, 
namely P = 5A3GeV/fm 3 and P 1 = 0.339GeV/ fm 3 . 
This corresponds to 2.495 x 10" 7 and 1.557 x 10" 8 in 
numerical units, respectively. The resulting pressure and 
velocity profiles for n/s = 0.01, are compared to the re- 
sults with previous LBS in Fig. [4] showing again very 
good agreement. It should be noted that, due to the 
above-mentioned pressure difference, the matter behind 
the shock moves with the velocity f3 ~ 0.6 (moderately 
relativistic regime), while the shock itself goes with the 
velocity v shock ~ 0.92. 

Note that the proposed model in the non-extended 
form (hereafter basic model) becomes numerically un- 
stable for higher velocities (/3 > 0.6). This instability 
is due to compressibility effects. It is known that the 
lattice Boltzmann method is intrinsically suited to low 
Mach number flows (low compressibility effects). There- 
fore, in order to overcome this problem, we use our ex- 
tended model, which enhances the stability of the nu- 
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x(fm) 



FIG. 4. Comparison between the current model and the pre- 
vious LBS model at rj/s = 0.01, for the pressure (top) and ve- 
locity (bottom) profiles in the moderately relativist ic regime. 



merical procedure without any appreciable loss of com- 
putational efficiency. To further investigate this issue, we 
carry out two simulations with the same conditions for 
relatively higher velocity, one using the basic model, and 
the other one using the extended one. The pressure is set 
to be P = 5A3GeV/fm 3 and P 1 = O.W95GeV/ fm 3 for 
the left side and the right side, respectively, which cor- 
responds to 2.495 x 10 _T and 7.785 x 10 -9 in numerical 
units. Here, rj/s = 0.01, S t /S x = 0.25, and a = 0.15 for 
the extended model. 

The results for the pressure and velocity profiles are 
shown in Fig. [5] at time t = 3.2 fm/ s. Note that the ap- 
plied pressure difference corresponds to f3 ~ 0.7. Using 
the basic model, an artificial discontinuity is observed 
in both the pressure and velocity profiles, which is due 
to instability problems in the numerical scheme. How- 
ever, the extended model proves capable of handling the 
simulation, the problem of the artificial discontinuity be- 
ing solved completely. Additionally, apart from the re- 
gion affected by the discontinuity, the good agreement 
between the results of two models can be interpreted as 
a validation for the precision of the extended model. The 
required CPU time for the simulation using the extended 
model and for the chosen value of S t /S x is 1069 ms. 

To the best of our knowledge, to date, there was no 
reported simulation of shock wave in viscous flow for 




1 2 3 4 5 6 

x (fm) 



FIG. 5. Comparison between the basic and extended models 
for the pressure (top) and velocity (bottom) profiles in the 
moderately relativistic regime. 



j3 > 0.6. However, for the inviscid case, which corre- 
sponds to the Euler equation at macroscopic scale, there 
exists an analytical solution for the Riemann problem. 
Therefore, in order to validate our extended model, we 
compare our results with the analytical solution of the in- 
viscid case in Ref. [35]. Hence, we need to solve the Euler 
equation by ignoring the viscous effects. It is worth men- 
tioning that, in the classical lattice Boltzmann method, 
the numerical solution becomes unstable as one tries to 
set the shear viscosity to zero (r = 1/2). This is also the 
case for our basic model. However, in the extended model 
we can solve the Euler equation by changing the collision 
step, such that instead of implementing the regular colli- 
sion, we simply set the discretized distribution functions 
to their corresponding equilibrium values. This is simi- 
lar to the procedure used in Ref. [36 to solve Euler equa- 
tion in the non-relativistic case. As one can notice from 
the Boltzmann equation, this corresponds to ignore the 
non-equilibrium part of the distribution function which 
contains the information about the dissipation. There- 
fore, we neglect the viscous effects in the macroscopic 
equations, obtaining the Euler equations. 

The results of our inviscid simulation are compared 
to the analytical results in Fig. [6] To drive the shock 
at velocity f3 ~ 0.9 (highly relativistic regime), the ap- 
plied pressure is set to Pq = 5A3GeV/ fm 3 and Pi = 
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Extended Model 
Analytical - 




Inviscid - 
r|/s=0.005 
r|/s=0.01 - 




FIG. 6. Comparison between the extended model and an- 
alytical results for the pressure (top) and velocity (bottom) 
profiles in the highly relativistic regime. 



FIG. 7. Comparison of the results for the inviscid case and 
different rj/s, for the pressure (top) and velocity (bottom) 
profiles in the highly relativistic regime. 



5A3MeV/fm 3 , which corresponds to 2.495 x 10 -7 and 
2.495 x 10 -10 in numerical units, respectively. The results 
are shown at t = 2.0 fm/s and very good agreement is 
found. The small discrepancies between our simulation 
and the analytical curve are related to the fact that a 
small value of bulk viscosity inevitably remains in our 
simulation, due to the coefficients, which are needed 
to increase the stability of the model. This issue makes a 
subject for future investigations. The CPU time required 
for the simulation shown in Fig. [6] is 593 ms. 

The fact that we can model properly the Euler equa- 
tion at very high velocities opens the possibility of using 
our model in astrophysical applications, where velocities 
are usually high and viscous effects are negligible. 

Fig. [7] shows the same result at f3 ~ 0.9 for different 
viscosities, compared to the inviscid case at t = 2.0 fm/s. 
The same conditions as mentioned above are considered 
again. The effects of increasing rj/s at this velocity are 
similar to the case of lower velocity. 

In order to demonstrate the ability of the model to 
simulate ultra-high velocities, we consider the case of 
the initial pressures Po = 5A3GeV/ fm 3 and Pi = 
0.0543M eF//ra 3 which corresponds to 2.495 x 10~ 7 and 
2.495 x 10 -12 in numerical units, respectively. We set 
6 t /S x = 0.25, a = 0.2, and rj/s = 0.01. The re- 
sults for the velocity profile and local Lorentz's factor 



at t = 2.0 fm/s are presented in Fig.[8j which shows that 
for this case j3 ~ 0.99 and j(u) ~ 9. This indicates that 
our model is numerically stable for simulating relativistic 
fluids with ultra- high velocities. The required CPU time 
for this simulation is 781 ms. 



Astrophysical application 

As an astrophysical application, we simulate a rela- 
tivistic shock wave, generated by, say, a supernova ex- 
plosion, colliding with a massive interstellar cloud, e.g. 
molecular gas [37]. The ejecta from the explosion of such 
supernovae are known to sweep the interstellar material 
up to relativistic velocities along the way (relativistic out- 
flows) [38]. We perform a three-dimensional simulation 
of a shock wave passing through a cold spherical cloud in 
a lattice of 200 x 100 x 100 cells. As mentioned earlier, 
in order to solve the equation of conservation of particle 
density, an extra distribution function is used, Eq.(58). 
As initial condition, the region is divided in two zones 
by the plane x = 65; at the left hand side (x < 65), 
the density is set to ri\ = 0.6cm -3 and the tempera- 
ture to Ti = 10 A K. The massive cloud is modeled as a 
solid sphere with radius of 10 cells and centered at the 
location (100,50,50), where we neglect the drag force 
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FIG. 8. Results of the simulation for the velocity profile (top) 
and Lorentz's factor (bottom) in the ultra-high relativistic 
regime. 



acting on the cloud due to the flow (the sphere will re- 
main at the same position during the whole simulation). 
Open boundary conditions are applied to all the external 
boundaries, except the left one, where an inlet boundary 
condition is applied by fixing the distribution function 
with the equilibrium distribution calculated at the initial 
condition, i.e. no and T . On the surface of the cloud, 
the cells evolve to the equilibrium distribution function 
evaluated at the constant values of n = ni, T = T\ and 
u = 0. It should be mentioned that, at each cell, the 
pressure can be calculated using the relation P = nT. 

By changing the initial condition at the right hand 
side of the dividing plane (x > 65), we are able to tune 
the velocity of the shock wave. Two different velocities 
are considered here. For the first case, we consider a 
shock wave at weakly relativistic regime (/? = 0.2), set- 
ting Tq = 2.71Ti and no = 0.9ni. For the second case, 
highly relativistic regime (/? = 0.9), theses values are 
T = 54.77Ti and n = 10.95ni. Fig.[9]shows the results 
of the 3D simulation of the shock wave, after colliding 
with the massive interstellar cloud for the density field 
for both cases. The simulations are performed in the in- 
viscid case, where a — 0.2 and S t /S x = 0.15. The density 
field is plotted in logarithmic scale at t = 1000 time step, 
where red and blue denote high and low values, respec- 
tively, and streamlines represent the velocity field. In 




FIG. 9. Snapshots of the three-dimensional simulation of a 
relativistic shock wave colliding with a massive interstellar 
cloud. Here, the density field is plotted in logarithmic scale 
in the weakly relativistic regime (top) and highly relativistic 
regime (bottom) at time t = 1000. The iso-surface in the 
second figure illustrates a region of low density (log(n/no) ~ 
-2.5). 



this figure, we observe an increase in the density down- 
stream of the collision, likely due to the fact that, during 
the propagation, the shock wave collects interstellar ma- 
terial and pushes it against the cloud (sweeping effect). 
However, at higher speeds, this increase on the density 
becomes less pronounced, and the passing of the shock 
wave generates a ring-shape region of low density down- 
stream (see isosurface in Fig. [9|. 

In order to study the viscous effects, the same sim- 
ulations are performed by introducing the dissipation 
and taking r = 1. Fig. [10| shows the pressure, density 
and temperature profiles at weakly and highly relativis- 
tic regimes for both, inviscid and viscous cases. The re- 
sults are shown along the x axis at y = z = 50 and at 
t = 600 time steps. Note that since Po, n o, an d Tq take 
very different values for the weakly and highly relativistic 
regimes, in order to make a clear comparison of the re- 
sults, we have normalized P, n, and T with Po, no, and 
Tq, respectively. As it can be appreciated, the viscous 
effects become relatively more important downstream of 
the cloud. Furthermore, it causes a pressure drop and 
decreases the density, while inducing a corresponding in- 
crement of temperature in the gas. Note, that the effects 
on the temperature are more significant at highly rela- 
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FIG. 10. Pressure (top), density (middle) and temperature 
(bottom) profiles in the weakly and highly relativist ic regimes, 
for the inviscid and viscous case. The results are shown along 
the x axis at y = z = 50 at t — 600 time steps. The thicker 
lines denote the regions where the interstellar cloud is located. 



tivistic regime than the pressure drop, while an opposite 
behavior is found at weakly relativist ic regime. 



IV. CONCLUSIONS 

In this paper, we have introduced a relativistic lat- 
tice Boltzmann model that is able to handle relativistic 
fluid dynamics at very high velocities. For this purpose, 



we have first expanded the Maxwell Jiittner distribution 
in orthogonal polynomials, by assuming as weight func- 
tion the equilibrium distribution at the local rest frame. 
A discretization procedure has been applied in order to 
adjust the expansion to the D3Q19 cell configuration, 
which, in order to avoid a multi-time evolution of the 
Boltzmann equation, leads to the problem of recovering 
only the conservation of the momentum-energy tensor. 
However, in the ultra-relativistic regime (e = Sp) the en- 
tire dynamics of the system is governed by this equation 
and the first order moment is not required. To extend 
the model to high velocities (/? ~ 1), we use a flux lim- 
iter scheme and introduce a bulk viscosity term into the 
Boltzmann equation, to increase the numerical stability 
in presence of discontinuities. 

In order to validate our model, we have compared the 
numerical results for shock waves in viscous quark-gluon 
plasmas with the results of other existing models, and 
found very good agreement. In addition, to the best of 
our knowledge, we have for the first time successfully sim- 
ulated shock waves in relativistic viscous flow for f3 > 0.6. 
We have also suggested a way to simulate near-inviscid 
flows (Euler equation) using the extended model by mod- 
ifying the collision step. For this case, we have compared 
the results with the analytical solution, finding again very 
satisafctory agreement. This offers a promising strategy 
to study astrophysical flows at very high speeds and neg- 
ligible viscous effects. Additionally, we have shown that 
our model is capable of simulating the Riemann prob- 
lem at ultra-high relativistic flows (7 ~ 10). Finally, we 
have studied the collision of a shock wave colliding with a 
massive interstellar cloud in weakly and highly relativis- 
tic regimes, for both inviscid and viscous cases. 

Summarizing, we have proposed a model which is 
simple, numerically efficient, and capable of simulating 
highly relativistic flows. It can also be used to simulate 
near-inviscid flows. Moreover, like other lattice Boltz- 
mann methods, our model is highly adaptable to paral- 
lel computing and can be used to simulate complex ge- 
ometries. These features make this model appealing for 
prospective applications in relativistic astrophysics and 
high energy physics. Extensions of this model to include 
higher order lattices so as to recover more moments of the 
equilibrium distribution, make a very interesting subject 
for future research. 
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